for i=1: NIut
    for j=1:NJut
        if (i<=NIue && j>NJue) continue; end
        
        indexU_cen   = Umap (i  ,j  );
        
        if (j<=NJue)
            
            % Boundary conditions on inlet
            if (i==1)
                Dp (indexU_cen,:) = 0;
                Dp (indexU_cen,indexU_cen) = 1;
                continue
            end
            
        else
            
            % Boundary conditions on expansion & outlet
            % Boundary condition expansion
            if (i==NIue+1)
                Dp (indexU_cen,:) = 0;
                Dp (indexU_cen,indexU_cen) = 1;
                continue
            end
            
        end
        % Boundary condition oulet - null normal gradient
        if (i==NIut)
            indexU_left  = Umap (i-1,j  );
            Dp (indexU_cen,:) = 0;
            Dp (indexU_cen,indexU_cen) = 1;
            Dp (indexU_cen,indexU_left) = -1;
            continue
        end
    end
end

for i=1:NIvt
    for j=1:NJvt
        if (i<=NIve && j>NJve) continue; end
        indexV_cen  = dimU + Vmap (i  ,j  ) ;
        
        if (j<=NJve)
            % Boundary conditions on top and bot boundaries
            if (j==1)
                Dp (indexV_cen,:) = 0;
                Dp (indexV_cen,indexV_cen) = 1;
                continue
            end
            if (j==NJve && i<=NIve )
                Dp (indexV_cen,:) = 0;
                Dp (indexV_cen,indexV_cen) = 1;
                continue
            end
        else
            % Boundary conditions on top and bot boundaries
            if (j==NJvt )
                Dp (indexV_cen,:) = 0;
                Dp (indexV_cen,indexV_cen) = 1;
                continue
            end
        end
        
    end
end